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CALCULATING R0T0RDYNAMIC COEFFICIENTS OF SEALS 
BY FINITE-DIFFERENCE TECHNIQUES 


F.J. Dletzen and R. Nordmann 
University of Kaiserslautern 
Kaiserslautern, Federal Republic of Germany 


For modelling the turbulent flow in a seal the Navier-Stokes equations in con- 
nection with a turbulence model (k-e -model) are solved by a finite-difference 
method. A motion of the shaft round the centered position is assumed. After cal- 
culating the corresponding flow field and the pressure distribution, the rotor- 
dynamic coefficients of the seal can be determined. These coefficients are com- 
pared with results obtained by using the bulk flow theory of Childs [1] and with 
experimental results. 


INTRODUCTION 

It is well known that the fluid forces in seals, which are described by equa- 
tion (1) 
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have a strong influence on the dynamic behaviour of rotating turbo-machinery. 

While there exist some good theories for calculating the coefficients of 
straight seals [1], no satisfactory model is known to describe the effects of 
grooved seals. Reference [2] presents a survey and comparison of results of 
existing theories. The authors' opinion is that the existing methods are not 
at all satisfactory. The main weakness of these theories is the fact, that they 
are using so called 'bulk-flow-theories' which connect the wall shear stress with 
the mean flow-velocity relative to this wall. However, in the region of a groove 
there occur stresses in the fluid which cannot be neglected. Calculating the flow 
by using the Navier-Stokes equations in connection with a turbulence model elimi- 
nates this disadvantage. Therefore, a finite difference model is presented which 
allows the calculation of the coefficients by using these equations. 
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Nomenclature: 


K, k 
D, d 
M, m 
U, V, w 

p 

k 

e 

u 6 ’ yp y t 
p 
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x, r, e 
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Dj,,, a » K 
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Forces on the shaft in z and y direction 
direct and cross-coupling stiffness in eq. (1, 24) 
direct and cross-coupling damping in eq. (1,24) 
direct and cross -coupl ing inertia in eq. (1, 24) 
axial, radial and circumferential velocity 
pressure 

turbulence energy 
energy dissipation 

effective, laminar and turbulent viscosity 

density 

time 

axial, radial and circumferential coordinate 
radial coordinate after transformation 
Constants od the k-e -model 
Constants of the k-£-model 

general variable standing for u, v, w, p, k or e 
general source term 

seal clearance by centric shaft position 
seal clearance by eccentric shaft position 
radius of the precession motion of the shaft 
perturbation parameter 
rotational frequency of the shaft 
precession frequency of the shaft 
entrance lost-coefficient 
Length of the seal 
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radius of the rotor (shaft) 
radius of the stator 


r i 


Subscripts 

0 

1 

R 

S 


zeroth order variables 
first order variables 
rotor 
stator 


MATHEMATICAL MODEL 


To describe turbulent flow by the Navier-Stokes equations the velocities and the 
pressure are separated into mean and fluctuating quantities. 

u=u+u' v=v+v' 

w=w+w' p=p+p' 


Time-averaging of the Navier-Stokes equations leads to terms of the following 
form: u'v 1 , v 'w' , u 'w' . 

To substitute these terms one can use the Boussinesq's eddy-viscosity concept. 
For example: 


u'v' 
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( 2 ) 


is the turbulent viscosity, which is not a fluid property but depends strongly 
on the state of flow. Summing up the laminar and turbulent viscosity to an effec- 
tive viscosity 

y e = y l + y t * (3) 

one obtains the following time-averaged Navier-Stokes equations for turbulent 
flow. (In the following the overbars are omitted.) 
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2. radial momentum: 
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3. tangential momentum 
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4. continuity equation 


|-x (pu > + rl'r (rpw) + rle (pw) = 0 


(7) 


To describe p t we use the k-e turbulence model [3, 4]. This model determines u t 
as a function of the kinetic energy k of the turbulent motion and the energy 
dissipation e. It is relative simple and often used to calculate the turbulent 
flow in seals [12, 13, 14, 15]. Stoff [12], for example, compares his flow meas- 
urements in a labyrinth seal with calculations on base of the k-e model. He 
observes that both agree well. 
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The equations for k and e can be derived in exact form from the Navier-Stokes 
equations 


5. turbulence energy k 
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6. energy dissipation e 
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To model the flow in the case of a shaft moving on an eccentric orbit, a coordi- 
nate-transformation [5, 6] is made. (Fig. 1) 


6 ( 0 , t) is the seal clearance, varying with angle 6 and time t. By this trans- 
formation the eccentric moving shaft is reduced to a shaft rotating in the 
centre of the seal. 

We must note that the following relations of the transformation must be used. 


PERTURBATION ANALYSIS 


The rotor dynamic coefficients of a seal are in a wide range independent of the 
shaft eccentricity. Therefore we can assume small shaft motions around the cen- 
tered position which allow us to use a perturbation analysis. 

6 = C 0 - ehj 

u = u + eu, v = v + ev, 

o 1 o I 

w = w Q + ew x P = P 0 + ep x 

With these expressions and the coordinate-transformation equation (12) the 
equations (4), (5), (6), (7), (9), (10) change themselves. 

This is demonstrated in the following examples. 

From equation (12) we obtain: 
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and so: 
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One obtains a set of zero-order equations for u Q , v Q , w Q , p Q , k Q , e Q and a set 
of first-order equations for Up Vp Wp pp It is assumed that the viscosity y g 
remains constant for small motions. Therefore the k^ and equations can be 
dropped. 

The variation of the seal clearance for an eccentric shaft can be described by 
the following equation 


h, = |cos0 + ^ sin0 
1 e e 


So we establish the same assumptions as in [1], that the velocities and the 
pressure in circumferential direction can be described by sin- and cos-functions, 
in our first order equations 

U 1 = U 1 C cos0 + u ls sin0 V 1 = v lc cos0 + v ls sin0 

W 1 = w lc cos0 + w ls sin0 P 1 = p lc cos0 + p ls sin0 

By separating in the resulting equations the terms with sin9 and cose we obtain 
two real equations of every 1. order equation. These equations are then arranged 
in a new form by introducing complex variables. 
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U 1 = u lc + iu ls 

*1 = w lc + iw ls 
hj_ = z + iy 


vi = v lc + iv ls 
Pi = Pic + ip ls 


We now assume that the shaft is moving on a circular orbit with frequency fi 
around the centered position. Also h^ takes the form: 



and similary 
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In the following we assume that t = 0; this means that the shaft is just moving 
through the z axis in the y-di recti on. 

The resulting equations for u Q , v 0 , w , p Q , k Q , e Q and Up Vp w-^ have all the 
same form. 
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Zeroth Order Equations 
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First Order Equations 

* r * 


ORIGINAL T 

OR POOR QUA 


rs rv" 


a i p e ~3x" 1+ 3x^ u e3x^ + n9n^ y e3x^ " 3x^ pu o 0 l^ " iW npu oV 


*1 p e 




b. 


w„ 


V°1 + ' p (- a + ir )0 l + 'n u o*l ' , S e H 1 + D 1 + i0 2 


-g- 1 * S*^) * fe<™«.g4 - fxlPV!) - &<■»»„»,) 

- ip eli<T> * $»1 + ip <^ ' °)*1 * < 2 n w o + i2 ^ + 1 n v o )fl l + °3 + 1D 4 

# - 4<r°i> ' ^ Vi> - l;< ow o Q i ) - Mi<"Vl> 

y G y a ft i a ny p 

- 13 ^*1 - <n v o + ip “ + ' 12 n"o> # l + D 5 + iD 6 


Only the first order continuity-equation to determine shows a slightly 
modified form. 


fj(pQl) - ^(nP«!) = 1®#! * t> 7 +fD 8 (20) 

The parameters Dj - Dg do not depend on Up Vp Wp p-^ and result from the 
coordinate-transformation. (D^ - Dg are shown in the appendix.) 

FINITE-DIFFERENCE METHOD 

For solving these equations a finite-difference procedure is used which is 
based on a method published by Gosman and Pun [7]. The seal is discretized 
by a grid (Fig. 2) and the variables are calculated at the nodes. The veloci- 
ties u, v are determined at points which lie between the nodes where the vari- 


85 


ables p, w, k, e are calculated (Fig. 3). Because of its general convergence 
a 'hybrid - difference' method is used, which means that the convective terms 
are calculated by a 'upwind 1 - or a 'central -difference' method as a function 
of flow-vel ocity and grid-distance. 

Because there is no explicit equation to calculate p we use the continuity 
equation. Starting by a guess for p, the momentum-equations are solved; with 
the resulting values u, v (Up Vp w^) the flow through the control-area around 
a point for the pressure p (p,) is calculated. If the difference between the 
entrance- and the exit-flow rate is less than 0, p must be reduced; in the 
opposite case p must be increased. This is done by the 'SIMPLE '-procedure [8] 
or better by the more modern version ' P ISO ' [9]. 

However one has to respect in the determination of p^ with these procedures 
that the equation for w^ has not the same form as for u^ and Vp Also we have 
to notice that u Q , v Q , p Q , w Q , k Q , e o are real-, while u p Vp w p p x are of 
complex type. The mesh to calculate u, v doesn't extend all the way to the boun- 
dary wall, and the component u, w is allowed to slip in accordance with the 
logarithmic law of the wall. 


LEAKAGE FLOW AND DYNAMIC COEFFICIENTS 
Leakage flow. Centered Position 

For centered shaft position the values u Q , v Q , p Q , yi Q , k Q e Q are determined. 
Boundary conditions: 

u oS = 0 v oS = 0 W 0 S = 0 

u oR = 0 v oR s 0 w oR = “- r i 

The leakage results from the calculated axial velocity u Q . 

Dynamic Coefficients, Eccentric Shaft Motion 

For calculating the dynamic coefficients the following assumptions are made: 

1. The shaft rotates on a circular orbit around the centered position. 

2. At time t = 0 the shaft is located at: z = r , v = 0 

0 J 

3. The viscosity remains constant in spite of the eccentric motion. 
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Boundary conditions: (Fig. 4, 5) 

Stator: u^ = (0., 0.) v^ = (0., 0.) = (0., 0.) 

Rotor: ^ 1 R = ^ 1 R ~ (0. ,(^-u))*Cq) = 0*) 

1 2 

Entrance: P/\ = 7 pu (1 + £) + P B 

Pl Bj = ' p u °Bj ^ 1 + ?) °l Bj 
Exit: P 1 q . = (0. , 0. ) 


To satisfy the entrance condition we make use of the iterative character of the 

finite-difference method. This means that we start with a pressure pig. at the 

entrance and after every iteration step we check if the calculated uig J satisfy 

condition (21). If not, the pressure pi R will be corrected. 

D j 


The resulting forces on the shaft are calculated by a pressure integration for 
the five precession frequencies: ft = 0u>, 0,5 00 , 1,0 00 , 1,5 w, 2 w. 


■ F z = V- 

/ pig dx 

( 22 ) 

TTr • 


/ pi^ dx 

(23) 


By a 'Least-Square-Fit' we obtain the rotor-dynamic coefficients of (1) from the 
following equations 

- F, = K + CS2 - MS 2 

z (24) 

- Fy = -k + Cft + mft^ 

The precession frequencies can be arbitrarily chosen, because the dynamic coeffi- 
cients are mostly independent of them. We take the same as in [1]. 


RESULTS FOR AN ANNULAR SEAL 

To test the theory, calculations are made for a straight smooth seal. The results 
are compared with the experimental values of Massmann [10] and the results of 
Childs theory [1]. 


Seal Data: 


L = 23,5 mm 

r- = 23,5 ram 

C = o,2 mm 
o 


U 1 = o,7-lo" 3 Ns/m 3 

p = 996 kg/m 3 

K -o,5 


A fully developed turbulent axial and circumferential flow at the entrance of the 
seal is assumed. As in [10] flowrates are measured, in the presented calcu- 
lation we suppose that the axial flow velocity is known: 


u average = 16,46 

and that the average circumferential velocity at the entrance is half the shaft- 
speed. 


For a known mass flow the pressure difference between entrance and exit of the 
seal can be calculated. The results of this theory are compared with Childs 
theory in Fig. 6. 

In Fig. 7, 8, 9, 10, 11 results of the presented theory, Childs theory [1] and 
experimental data from Massmann [10] are shown. 


Both theories are in good agreement with each other and with the measurements. 
For calculation a mesh with 15 x 5 nodes in x-r direction was applied. The CPU 
time was about 30 sec on a Siemens 7.561 computer. 


RESULTS FOR A GROOVED SEAL 

We also made some calculations, for the grooved seal, whose geometry and seal 
data are shown in Fig. 12. In Fig. 13 the leakage for a given pressure difference 
is presented as a function of the groove depth. First the leakage decreases and 
then slightly increases again. This behaviour agrees with the measurements of 
Black [11]. 

In Fig. 14, 15, 16, 17 the stiffness K, k and the damping D, d are shown. The 
coefficients K, k, D decrease with growing groove depth. Only the damping d 
increases. Although we haven't yet any experimental results for this seal the 
tendencies seem to be right. 


CONCLUSION 

It is shown that it is possible to calculate the dynamic coefficients of seals 
with a finite-difference method, based on the Navier-Stokes equation in connec- 
tion with a turbulence model. Although application on straight seals is possible, 
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it was not our aim to develop a procedure for this seal configuration but to 
present a method which will be applicable on grooved seals. The superiority of 
the theory versus other methods is the simplicity in use for grooved seals by 
only neglecting the terms - Dg in the equations for Dp Vp Wp Pj in the 
grooves, while there exists no mesh displacement. 

Appendix: 

Transformation-constants for first order equations 
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Fig. 2 Mesh arrangement in the seal 
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Fig. 6 Pressure difference as a function of the shaft 
speed for an annular seal 
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Fig. 7 Direct-inertia as a function of the shaft speed 
for an annular seal 
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Fig. 8 Direct damping as a function of the shaft speed 
for an annular seal 
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Fig. 9 Cross-coupling damping as a function of the 
shaft speed for an annular seal 
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Fig. 14 Direct stiffness of grooved seals as a Function 
of the groove depth 
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Fig. 15 Cross-coupling stiffness of grooved seals as a function 
of the groove depth 
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Fig. 16 Direct damping of grooved seals as a function of the groove depth 


d 

(Ns/m) 

80 


60 


40 


20 


0 0.5 1.0 1.5 2.0 H r /C q 



Fig 17 Cross-coupling damping of grooved seals as a function 
of the groove depth 
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